#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>
#include <cufft.h>
#include "prof.h"
#include <pthread.h>
#include <unistd.h>
// CUDA错误检查宏
#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            fprintf(stderr, "CUDA错误 在 %s:%d 错误码=%d(%s) 调用: %s\n", \
                    __FILE__, __LINE__, err, cudaGetErrorString(err), #call); \
            exit(EXIT_FAILURE); \
        } \
    } while (0)

// cuFFT错误检查宏
#define CUFFT_CHECK(call) \
    do { \
        cufftResult err = call; \
        if (err != CUFFT_SUCCESS) { \
            fprintf(stderr, "cuFFT错误 在 %s:%d 错误码=%d\n", \
                    __FILE__, __LINE__, err); \
            exit(EXIT_FAILURE); \
        } \
    } while (0)

// 上采样内核：将输入信号x_d上采样p倍，生成x_up_d
__global__ void upsample_kernel(float *x_up_d, float *x_d, int N, int p) {
    int m = blockIdx.x * blockDim.x + threadIdx.x;
    int L = N * p; // 上采样后信号长度
    if (m < L) {
        if (m % p == 0) { // m是p的倍数时，复制原始样本
            int n = m / p;
            x_up_d[m] = x_d[n];
        } else { // 其余位置填充0
            x_up_d[m] = 0.0f;
        }
    }
}

// 滤波内核：在频域中应用低通滤波器
__global__ void filter_kernel(cufftComplex *Y_d, cufftComplex *X_up_d, int cutoff, int L) {
    int k = blockIdx.x * blockDim.x + threadIdx.x;
    if (k <= L/2) { // 实数FFT输出长度为L/2 + 1
        if (k <= cutoff) { // 保留截止频率内的频谱
            Y_d[k] = X_up_d[k];
        } else { // 高于截止频率的部分置零
            Y_d[k].x = 0.0f;
            Y_d[k].y = 0.0f;
        }
    }
}

// 缩放内核：cuFFT的IFFT结果需要除以信号长度L
__global__ void scale_kernel(float *y_up_d, int L) {
    int m = blockIdx.x * blockDim.x + threadIdx.x;
    if (m < L) {
        y_up_d[m] /= (float)L;
    }
}

// 下采样内核：从y_up_d中每隔q个样本取一个，生成y_d
__global__ void downsample_kernel(float *y_d, float *y_up_d, int q, int L, int M) {
    int k = blockIdx.x * blockDim.x + threadIdx.x;
    if (k < M) {
        int idx = k * q;
        if (idx < L) { // 确保索引不超过上采样信号长度
            y_d[k] = y_up_d[idx];
        } else {
            y_d[k] = 0.0f; // 超出部分填充0
        }
    }
}

int main(int argc, char** argv) {
    // --- 参数设置 ---
    int N = 1024000; // 输入信号长度
    int p = 2;    // 上采样因子
    int q = 3;    // 下采样因子
    int L = N * p; // 上采样后信号长度
    int M = (N * p + q - 1) / q; // 输出信号长度，ceil(N*p/q)

    // 获取运行次数n
    if (argc != 2) {
        fprintf(stderr, "用法: %s <运行次数>\n", argv[0]);
        return 1;
    }
    int n = atoi(argv[1]);
    if (n <= 0) {
        fprintf(stderr, "运行次数必须为正整数。\n");
        return 1;
    }

    profInit();//初始化
    profEnterSect(0);
    profSetSectName(1, "resamole htod_time ");//根据需要设置计时段名称
    profSetSectName(2, "resamole compute_time");
    profSetSectName(3, "resamole dtoh_time");

    // --- 分配主机内存 ---
    float *x = (float*)malloc(sizeof(float) * N); // 输入信号
    float *y = (float*)malloc(sizeof(float) * M); // 输出信号

    // 生成随机输入数据
    for (int i = 0; i < N; i++) {
        x[i] = (float)rand() / RAND_MAX; // 随机数范围[0,1]
    }

    // --- 分配设备内存 ---
    float *x_d, *x_up_d, *y_up_d, *y_d; // 输入、上采样、中间结果、输出
    cufftComplex *X_up_d, *Y_d; // FFT和滤波后的频谱
    CUDA_CHECK(cudaMalloc(&x_d, sizeof(float) * N));
    CUDA_CHECK(cudaMalloc(&x_up_d, sizeof(float) * L));
    CUDA_CHECK(cudaMalloc(&X_up_d, sizeof(cufftComplex) * (L/2 + 1)));
    CUDA_CHECK(cudaMalloc(&Y_d, sizeof(cufftComplex) * (L/2 + 1)));
    CUDA_CHECK(cudaMalloc(&y_up_d, sizeof(float) * L));
    CUDA_CHECK(cudaMalloc(&y_d, sizeof(float) * M));

    // --- 创建cuFFT计划 ---
    cufftHandle plan_r2c; // 实数到复数FFT计划
    CUFFT_CHECK(cufftPlan1d(&plan_r2c, L, CUFFT_R2C, 1));
    cufftHandle plan_c2r; // 复数到实数IFFT计划
    CUFFT_CHECK(cufftPlan1d(&plan_c2r, L, CUFFT_C2R, 1));

    // --- 计算滤波器截止频率 ---
    double min_pq = fmin(1.0/p, 1.0/q); // min(1/p, 1/q)
    int cutoff = (int)floor(min_pq * L / 2.0); // 截止频率对应的频点索引

    // --- 设置线程块和网格大小 ---
    int block_size = 256;
    int grid_size_upsample = (L + block_size - 1) / block_size;
    int grid_size_filter = ((L/2 + 1) + block_size - 1) / block_size;
    int grid_size_downsample = (M + block_size - 1) / block_size;

    // --- 预热运行（不计入统计） ---
    CUDA_CHECK(cudaMemcpy(x_d, x, sizeof(float) * N, cudaMemcpyHostToDevice));
    upsample_kernel<<<grid_size_upsample, block_size>>>(x_up_d, x_d, N, p);
    CUFFT_CHECK(cufftExecR2C(plan_r2c, x_up_d, X_up_d));
    filter_kernel<<<grid_size_filter, block_size>>>(Y_d, X_up_d, cutoff, L);
    CUFFT_CHECK(cufftExecC2R(plan_c2r, Y_d, y_up_d));
    scale_kernel<<<grid_size_upsample, block_size>>>(y_up_d, L);
    downsample_kernel<<<grid_size_downsample, block_size>>>(y_d, y_up_d, q, L, M);
    CUDA_CHECK(cudaMemcpy(y, y_d, sizeof(float) * M, cudaMemcpyDeviceToHost));

    // --- 执行n次运行并计时 ---
    for (int i = 0; i < n; i++) {
        // 记录htod时间
        profEnterSect(1);
        CUDA_CHECK(cudaMemcpy(x_d, x, sizeof(float) * N, cudaMemcpyHostToDevice));
        CUDA_CHECK(cudaDeviceSynchronize());
        profLeaveSect(1);

        // 记录计算时间
        profEnterSect(2);
        upsample_kernel<<<grid_size_upsample, block_size>>>(x_up_d, x_d, N, p);
        CUFFT_CHECK(cufftExecR2C(plan_r2c, x_up_d, X_up_d));
        filter_kernel<<<grid_size_filter, block_size>>>(Y_d, X_up_d, cutoff, L);
        CUFFT_CHECK(cufftExecC2R(plan_c2r, Y_d, y_up_d));
        scale_kernel<<<grid_size_upsample, block_size>>>(y_up_d, L);
        downsample_kernel<<<grid_size_downsample, block_size>>>(y_d, y_up_d, q, L, M);
        CUDA_CHECK(cudaDeviceSynchronize());
        profLeaveSect(2);

        // 记录dtoh时间
        profEnterSect(3);
        CUDA_CHECK(cudaMemcpy(y, y_d, sizeof(float) * M, cudaMemcpyDeviceToHost));
        CUDA_CHECK(cudaDeviceSynchronize());
        profLeaveSect(3);

        // 同步并计算时间
        // CUDA_CHECK(cudaEventSynchronize(stop_dtoh));

    }




    profLeaveSect(0);
    profPrintInfo();

    // --- 释放资源 ---
    CUDA_CHECK(cudaFree(x_d));
    CUDA_CHECK(cudaFree(x_up_d));
    CUDA_CHECK(cudaFree(X_up_d));
    CUDA_CHECK(cudaFree(Y_d));
    CUDA_CHECK(cudaFree(y_up_d));
    CUDA_CHECK(cudaFree(y_d));
    CUFFT_CHECK(cufftDestroy(plan_r2c));
    CUFFT_CHECK(cufftDestroy(plan_c2r));
    free(x);
    free(y);



    return 0;
}

